Integrated omics reveal novel functions and underlying mechanisms of the receptor kinase FERONIA in Arabidopsis thaliana

Abstract The receptor kinase FERONIA (FER) is a versatile regulator of plant growth and development, biotic and abiotic stress responses, and reproduction. To gain new insights into the molecular interplay of these processes and to identify new FER functions, we carried out quantitative transcriptome, proteome, and phosphoproteome profiling of Arabidopsis (Arabidopsis thaliana) wild-type and fer-4 loss-of-function mutant plants. Gene ontology terms for phytohormone signaling, abiotic stress, and biotic stress were significantly enriched among differentially expressed transcripts, differentially abundant proteins, and/or misphosphorylated proteins, in agreement with the known roles for FER in these processes. Analysis of multiomics data and subsequent experimental evidence revealed previously unknown functions for FER in endoplasmic reticulum (ER) body formation and glucosinolate biosynthesis. FER functions through the transcription factor NAI1 to mediate ER body formation. FER also negatively regulates indole glucosinolate biosynthesis, partially through NAI1. Furthermore, we found that a group of abscisic acid (ABA)-induced transcription factors is hypophosphorylated in the fer-4 mutant and demonstrated that FER acts through the transcription factor ABA INSENSITIVE5 (ABI5) to negatively regulate the ABA response during cotyledon greening. Our integrated omics study, therefore, reveals novel functions for FER and provides new insights into the underlying mechanisms of FER function.

FER plays critical roles in plant growth and stress responses. The functional and mechanistic studies of FER can illustrate how a cell surface receptor kinase integrates diverse signals to profoundly influence various cellular pathways and thus balance growth and stress in response to environmental conditions. Studies on FER also provide an important paradigm for a better understanding of receptor kinase signaling in general.
To gain new insights into the functions and mechanisms associated with FER, we globally quantified transcripts, proteins, and phosphorylation sites in wild-type (WT) plants and the fer-4 mutant. Our study validated many known functions of FER and revealed several interesting novel functions for this receptor kinase. We discovered that FER is involved in endoplasmic reticulum (ER) body formation and glucosinolate biosynthesis. FER functions through the transcription factor NAI1 to negatively mediate ER body formation. FER also negatively regulates indole glucosinolate biosynthesis, partially through NAI1. We also determined that FER phosphorylates and destabilizes abscisic acid (ABA)-INSENSITIVE 5 (ABI5), a critical transcriptional regulator of cotyledon greening (Guan et al., 2014), which provides a novel mechanism underlying the negative role of FER in ABA signaling. Taking advantage of the comprehensive integrated omics analysis along with genetics, cell biology, and biochemistry experiments, our study reveals novel functions for FER and provides new insights into the underlying mechanisms of FER function.

Integrated omics confirm known functions and predict novel functions of FER
We performed multiomics analyses on 4.5-week-old Arabidopsis rosette leaves from the WT Columbia-0 (Col-0) and fer-4 (hereafter fer) grown under long-day conditions (see "Materials and methods"). Under our growth conditions, 4.5-week-old fer plants showed a pronounced stunted growth phenotype with smaller rosettes, while 10-day-old fer seedlings exhibited an overall similar growth to that of WT (Supplemental Figure S2, C and D). Transcriptome profiling by 3 0 -based RNA sequencing (QuantSeq) identified 3,908 (q 5 0.05) transcripts that are differentially expressed in fer (Supplemental Figure S1), of which 2,299 had increased levels and 1,609 decreased levels in the mutant compared to WT ( Figure 1A; Supplemental Data Set 1). More than 50% of the 3,908 genes (2,190/3,908) were also identified as differentially expressed in fer mutants in our previous report (Guo et al. 2018;Supplemental Figure S2B). We determined that 4,129 (q 5 0.05) proteins have different abundance in the fer mutant out of 8,621 quantified proteins (i.e. protein groups; Supplemental Figure S1). Among the differentially abundant proteins (DAPs), 2,346 had increased and 1,783 had decreased levels in the fer mutant ( Figure 1A; Supplemental Data Set 1). Finally, we detected 11,955 phosphosites from 2,959 proteins (i.e. protein groups), among which 3,432 phosphosites (q 5 0.2) had altered phosphorylation status in the mutant background relative to the WT (Supplemental Figure S1). Further, 1,577 phosphosites (from 703 different proteins) had elevated phosphorylation, and 1,855 phosphosites (from 691 different proteins) exhibited decreased phosphorylation, with 110 proteins displaying both elevated and decreased phosphorylation levels at different sites ( Figure 1A; Supplemental Data Set 1). Comparisons among the three sets of omics data revealed that 29% of all DAPs also exhibit differential expression of their cognate transcripts in fer, with most of the DAPs following the same direction of change as their transcript ( Figure 1B; Supplemental Figure S2A). Similarly, 47% of the differentially phosphorylated proteins also showed an alteration at the corresponding transcript or protein levels ( Figure 1B; Supplemental Figure S2A).
We performed gene ontology (GO) analysis using ClueGO in Cytoscape on our multiomics data (Bindea et al., 2009) (Supplemental Figures S3-S8;Supplemental Data Set 2). We considered GO terms with a corrected P 5 0.05 as significantly enriched. The GO analysis validated many previous findings concerning FER function in plants. Namely, GO terms related to stress phytohormones were enriched ( Figure 1C; Supplemental Data Set 2). For example, GO terms such as response to JA (GO:0009753), response to ABA (GO:0009737), and ethylene (GO:0071369), were significantly enriched among transcripts and/or proteins with increased levels in fer, which corroborated the previous findings that FER regulates ABA, JA, and ethylene signaling pathways and that loss-of-function fer mutants are hypersensitive to these phytohormones (Deslauriers and Larsen, 2010;Yu et al., 2012;Guo et al., 2018). Additionally, GO terms related to growth-promoting phytohormones were enriched ( Figure 1C). Consistent with a role for FER in brassinosteroid (BR)-mediated plant growth, the GO terms response to BR (GO:0009741) and BR-mediated signaling pathway (GO:0009742) were enriched among proteins with decreased phosphorylation levels (Guo et al., 2009b). Providing support for the findings that FER is involved in auxin-regulated processes, GO terms such as auxin polar transport (GO:0009926), auxin activated signaling pathway (GO:0009734), and response to auxin (GO:0009733) were also enriched among transcripts with lower levels in fer (Duan et al., 2010; Figure 1C; Supplemental Figures S3-S8; Supplemental Data Set 2).
We also observed the enrichment of GO terms related to various biotic stresses, for example, defense response to bacterium (GO:0042742), and defense response to fungus (GO:0050832; Figure 1C; Supplemental Data Set 2), which confirmed the involvement of FER in fungal and bacterial defense responses (Masachis et al., 2016;Stegmann et al., 2017;Guo et al., 2018). Similarly, GO terms related to abiotic stresses were enriched, for example, response to cold (GO:0009409), response to salt stress (GO:0009651), and response to osmotic stress (GO:0006970; Figure 1C; Supplemental Data Set 2). Consistent with our findings, FER has been experimentally shown to regulate responses to cold stress, salt stress, and osmotic stress imposed by mannitol. In agreement, the loss-of-function fer mutant is hypersensitive to cold and high salt and is resistant to osmotic stress Feng et al., 2018;Zhao et al., 2018).
Moreover, the analysis of the omics data predicted previously unknown functions for FER ( Figure 1C; Supplemental Data Set 2). Among these, ER body (GO:0010168) and glucosinolate metabolism (GO:0019760) were enriched, suggesting that FER is involved in ER body formation and glucosinolate metabolism.
Using the hypophosphorylated proteins in the fer mutant that are potential FER substrates, we predicted the consensus sequences of the FER phosphorylation sites using the motifeR software (Wang et al., 2019a), which identified 33 significantly enriched consensus sequences of FER phosphorylation motifs (Supplemental Data Set 3). The LOGOs of the three most enriched sequences are shown in Supplemental Figure S9A. Furthermore, we carried out an in vitro kinase assay to validate potential FER substrates. We selected nine proteins for this assay, each harboring at least one of the 33 enriched FER phosphorylation motifs. These proteins have diverse functions, including transcription factors (ABA-RESPONSIVE ELEMENTS-BINDING FACTOR3 [ABF3], OXIDATIVE STRESS2 [OXS2], FLOWERING BHLH3 [FBH3], ABI5, and GT2), a plasma membrane-localized protein (SHOU4) involved in cellulose biosynthesis (Polko et al., 2018), a cytoplasmic protein (NAIP1, At3g51950) involved in ER body formation (Wang et al., 2019b) and a nucleus/chloroplast-localized protein (NUCLEAR SHUTTLE INTERACTING, NSI). The identified phosphopeptides and phosphorylation sites in these proteins are shown in Figure 2A. We performed the above in vitro kinase assays with FER kinase (FERK, consisting of the FER C-terminus containing the kinase domain) fused to glutathione S-transferase (GST) as described (Guo et al., 2018), with the nine candidate substrate proteins fused to maltose-binding protein (MBP). The kinase assays revealed that seven of the nine proteins are phosphorylated in the presence of FERK ( Figure 2B), suggesting that FER can directly phosphorylate these proteins.
substrates, mFERK did not show appreciable levels of autophosphorylation or phosphorylation of any of these substrates when provided in amounts equal to those used for FERK ( Figure 2C). Interestingly, a larger amount of mFERK (i.e. approximately 200 times more than FERK) showed autophosphorylation in an in vitro kinase assay ( Figure 2D), suggesting residual kinase activity in this mutant. We estimated the amount of FERK and mFERK used for kinase assay by immunoblotting using an anti-FER antibody (Supplemental Figure S9B; Guo et al., 2018). It is worth Figure 2 FER directly phosphorylates many proteins with diverse functions. A, List of proteins selected for in vitro kinase assays, with gene identifier, phosphopeptide(s), phosphosites, and the change in phosphorylation levels in fer, as indicated by Log 2 (fold-change[fer/WT]). B, In vitro kinase assay of the selected proteins. Top, autoradiograph of the in vitro kinase assay, showing FERK autophosphorylation (red arrow) and phosphorylation of the selected proteins (asterisks). Note that NAIP1 phosphorylation overlaps with FER autophosphorylation. Bottom, SYPRO RUBY staining of the gel with the same amount of proteins run separately. C, In vitro kinase assay of FBH3, GT2, and NSI using a similar amount of FERK and FERK K565R (mFERK). D, Autophosphorylation of FERK and different amounts of mFERK. 1Â mFERK indicates the same amount as FERK. Note the change in mobility of FERK K565R compared to FERK when resolved on 8% SDS-PAGE. E, Protein-protein interaction by BiFC in N. benthamiana leaves. Fluorescence and light images of leaf epidermal cells are shown. Scale bars = 100 mm. F-G, Co-IP assays confirming the interactions of FER with SHOU4 (F) and FBH3 (G). Constructs were co-infiltrated in N. benthamiana. Proteins were immunoprecipitated with anti-FLAG (mouse) and detected with anti-FLAG (rabbit) and anti-GFP (rabbit) antibodies.
noting that the mFERK appears to migrate faster than FERK when resolved on 8% (w/v) sodium dodecyl-sulfate polyacrylamide gel electrophoresis (SDS-PAGE) (Supplemental Figure S9B). We also tested the potential for interaction between FER and its putative substrates using bimolecular fluorescence complementation (BiFC) assays in transiently infiltrated Nicotiana benthamiana leaves. FER interacted with SHOU4, FBH3, and ABI5 ( Figure 2E; Supplemental Figure S9C). DAWDLE (DDL), a mostly nucleus-localized protein that is not among the proteins detected in the phosphoproteomics data, was used as a negative control (Supplemental Figure  S9D). We confirmed the interaction between FER and the substrates SHOU4 and FBH3 by co-immunoprecipitation (Co-IP) assays in N. benthamiana leaves transiently infiltrated with constructs encoding tagged versions of each protein (Figure 2, F and G). The Co-IP did not detect an interaction between FER and ABI5, which might be due to a highly transient interaction and/or highly unstable FERphosphorylated ABI5.
Taken together, our multiomics data analysis validated previous findings that FER plays important roles in plant growth and development, abiotic and biotic stress responses. The phosphoproteomics detected many FER substrates in planta. We further validated some of the FER substrates using in vitro kinase assays, BiFC, and Co-IP. The multiomics data analysis also predicted novel functions for FER in ER body formation and glucosinolate biosynthesis, as well as FER regulation of ABI5 in ABA signaling. As detailed in the next sections, we explored and experimentally validated the role of FER in these processes.

FER negatively regulates ER body formation
The ER body is a type of membranous structure, 1 Â 10 mm in size, that is contiguous with the ER and is surrounded by ribosomes in the Arabidopsis cytoplasm (Hayashi et al., 2001). ER bodies, specific to Brassicales, are constitutively present in epidermal cells of healthy young seedlings but absent in rosette leaves of healthy WT plants (Hayashi et al., 2001). NAI1, a basic helix-loop-helix (bHLH) transcription factor, is required for ER body formation (Matsushima et al., 2003a(Matsushima et al., , 2004) and a loss-of-function nai1 mutant lacks ER bodies even in young seedlings. Many ER body-localized or body-related proteins are known (Matsushima et al., 2003b), and the expression of many of the genes encoding these proteins is regulated by NAI1 (Yamada et al., 2008).
Consistent with the enrichment of the GO term for ER body (GO:0010168) ( Figure 1C; Supplemental Data Set 2), 16 out of the 17 known ER body-associated genes (listed in Figure 3A; Wang et al., 2017) displayed increased transcript levels, and 14 of the16 had increased protein abundance in the fer mutant background compared to WT ( Figure 3A), suggesting that FER negatively regulates ER body formation. Sixteen of the 17 genes also exhibited increased transcript levels in the fer mutant in a previously described transcriptome dataset (fer-DEs-CB in Figure 3A) (Guo et al., 2018). We also generated an additional set of QuantSeq transcriptome data from 10-day-old seedlings for WT and the fer mutant and identified 7,718 differentially expressed transcripts (2,718 of them downregulated and 5,018 upregulated, respectively, in fer relative to WT) (Supplemental Data Set 4). All 17 ER body genes showed increased transcript levels ( Figure 3A). A volcano plot also revealed the significant alterations in abundance for their encoded proteins ( Figure 3B). We confirmed the expression pattern of three selected ER body genes by reverse transcription-quantitative polymerase chain reaction (RT-qPCR) in 3-week-old plants ( Figure 3E). Consistent with the transcriptome data from 4.5-week-old plants and 10-day-old seedlings, NAI1, NAI2, and PYK10 transcripts accumulated to higher levels in fer relative to WT. NAI2 and PYK10 were expressed at lower levels in the nai1 mutant, corroborating the previous findings that NAI1 activates NAI2 and PYK10 transcription (Yamada et al., 2008). To better understand the genetic interaction between FER and NAI1, we constructed the fer nai1 double mutant. Interestingly, the expression of NAI2 and PYK10 decreased in the fer nai1 double mutant compared to the fer single mutant, suggesting that FER regulates the expression of ER body genes through NAI1. It is worth noting that the nai1 mutant (CS69075) we used in this study has a point mutation (Matsushima et al., 2004) that does not disrupt the transcription of the mutant allele or its regulation by FER ( Figure 3E).
To visualize ER bodies in plants, we crossed the ER marker GFP-HDEL (encoding the green fluorescent protein [GFP] with the ER retention signal HDEL) into the fer mutant background. We observed ER bodies in 5-and 10-day-old seedlings and 3-week-old rosette leaves of fer GFP-HDEL (Figure 3, C-D; Supplemental Figure S10, A-C). We used 5day-old seedlings, in which ER bodies are abundant, and 10day-old seedlings, in which ER bodies are largely diminished in WT cotyledons expressing GFP-HDEL, in the rest of this study. In 5-day-old seedlings, we observed ER bodies in the cotyledons, hypocotyls, and roots of both GFP-HDEL and fer GFP-HDEL, while we detected no ER bodies in nai1 GFP-HDEL except for a very few in roots (Supplemental Figure  S10, B and C). We counted more ER bodies in the cotyledons and roots of fer GFP-HDEL compared to the GFP-HDEL line, likely due to increased NAI1 expression in the fer mutant (Supplemental Figure S10, B and C). While we failed to observe ER bodies in the cotyledons of 10-day-old GFP-HDEL seedlings ( Figure S10A). NAI1 was more highly expressed in fer relative to WT (Figure 3, A and E). To test our hypothesis that FER regulates ER body by negatively regulating NAI1, we generated the fer nai1 double mutant harboring the GFP-HDEL reporter, yielding fer nai1 GFP-HDEL. Similar to nai1 GFP-HDEL, 5-day-old fer nai1 GFP-HDEL seedlings lacked ER bodies (Supplemental Figure S10, B and C). The increased ER body  Figure 3 FER negatively regulates ER body formation. A, List of ER body-associated genes  and their regulation in the differentially expressed transcripts and DAPs in 4.5-week-old fer plants, differentially expressed transcripts in 10-day-old fer seedlings (Transcripts 10 days) and previously published fer transcriptome data (fer DEs (2018)) (Guo et al., 2018). Up: increased levels in fer mutant. NS: detected, but no statistically significant differences. ND: not detected. B, Volcano plot constructed with the 8,621 detected proteins, with ggplot2 (Wickham, 2016) in R using -Log 10 -transformed q-values and Log 2 (Fold-change) of protein abundance. The x-axis represents Log 2 (Fold-change) of protein levels between fer and WT; the y-axis represents statistical significance using -Log 10 (adjusted P-value). Black, no statistically significant changes in fer; blue, decreased levels in fer; red, increased levels in fer. Selected ER body-associated proteins are indicated. C, Number of ER bodies in 10-day-old cotyledons, quantified from confocal images of the four genotypes overexpressing the ER marker GFP-HDEL. WT (n = 31), fer (n = 18), nai1 (n = 31), and fer nai1 (n = 28). Data are shown as violin plots with median, first and third quartiles. Different letters indicate significant differences according to one-way analysis of variance (ANOVA)Tukey's multiple range tests (P 5 0.05). D, Representative confocal microscopy images from (C). Scale bar = 20 mm. E, RT-qPCR of representative ER body-related genes (NAI1, NAI2, and PYK10) in 3-week-old rosette leaves of WT, fer, nai1, and fer nai1 overexpressing GFP-HDEL. F, RT-qPCR of NAI1, NAI2, and PYK10 in 3-week-old WT, fer, fer FERpro:FER-GFP, or fer FERpro:FER K565R -GFP complementation lines. Relative expression levels were normalized to those of the reference gene ACTIN2. qPCR was performed with two to three technical replicates of three to four independent biological replicates. Data represent means ± SEM. Different letters indicate significant differences according to one-way ANOVA Tukey's multiple range tests (P 5 0.05).
phenotype observed in 10-day-old fer GFP-HDEL also diminished in the fer nai1 GFP-HDEL background ( Figure 3, C and D), suggesting that FER functions through NAI1 to regulate ER body formation. To further support a role for FER in ER body formation, we transformed an artificial microRNA (amiRNA) specific for FER (FER-amiRNA) into the GFP-HDEL reporter line to generate the FER-amiRNA GFP-HDEL line. The FER-amiRNA has previously been shown to effectively knock down endogenous FER transcripts, with the knockdown plants behaving similarly to fer mutants for both plant growth and bacterial defense in T 2 transgenic lines (Guo et al., 2009b(Guo et al., , 2018. In FER-amiRNA GFP-HDEL T 2 lines, FER protein abundance was much lower than in the GFP-HDEL line (Supplemental Figure S11A). We counted more ER bodies in 10-day-old cotyledons in all three FER-amiRNA GFP-HDEL T 2 lines tested here (lines #2, #4, and #8; Supplemental Figure S11, D and E). FER-amiRNA GFP-HDEL lines also exhibited a slight increase in the number of ER bodies in 5-day-old cotyledons compared to the GFP-HDEL line (Supplemental Figure S11, B and C), similar to fer GFP-HDEL (Supplemental Figure S10, B and C). In summary, these results suggest that FER negatively regulates ER body formation by negatively regulating NAI1.
Complementation assays with full-length intact FER and mutant FER K565R constructs revealed that FER kinase activity is not required for some FER-mediated processes, such as ovule fertilization (Kessler et al., 2015), but was for other processes, such as stomatal movement, vegetative growth, and root elongation Haruta et al., 2018). Chakravorty et al. further showed that high levels of FER K565R protein can complement the fer dwarf growth phenotype and restore stomatal movements regulated by the ligand RALF1. We thus investigated if ER body regulation by FER requires its kinase activity. From our transcriptome data and RT-qPCR results, we observed that the expression pattern of ER body genes is a good indicator of ER body phenotype. We measured the expression of ER body genes in WT, fer, and complementation lines harboring an intact FER or FER K565R transgene (fer FERpro:FER-GFP, fer FERpro:FER K565R -GFP #5 and fer FERpro:FER K565R -GFP #8; Shih et al., 2014;Chakravorty et al., 2018). Consistent with their reported growth phenotypes, the FERpro:FER transgene fully complemented the ER body gene expression pattern seen in fer. The lines FERpro:FER K565R -GFP #5 and #8 showed partial restoration, with line #8, with higher levels of the FER K565R -GFP fusion protein, showing greater rescue than line #5 ( Figure 3F). These results suggest that FER regulation of ER body formation is dependent on its kinase activity.
NAI2 is an ER body protein found only in Brassicales that is involved in ER body formation and function; loss-offunction nai2 mutants largely lack ER bodies (Yamada et al., 2008;Wang et al., 2019b). NAI1, the major transcription factor in ER body formation, was shown to activate NAI2 gene expression by directly binding to its promoter (Sarkar et al., 2021). To further support the conclusion that FER functions through NAI1 to regulate ER body formation, we generated FER-amiRNA nai2-2 GFP-HDEL lines, in which FER levels dropped relative to nai2-2 ( Figure 4A). Similar to nai2-2 GFP-HDEL, all three individual FER-amiRNA nai2-2 GFP-HDEL lines lacked ER bodies in the cotyledons of both 5-and 10-day-old seedlings, while FER-amiRNA GFP-HDEL showed ER body formation in both 5-and 10-day-old seedlings ( Figure 4, B-E). These results further demonstrate that FER negatively regulates NAI1 and hence NAI2 to control ER body formation.

FER negatively regulates indole glucosinolate biosynthesis
Glucosinolates are a family of secondary metabolites specific to Brassicales order that play important roles in responses to biotic stresses such as insect herbivory (Burrow et al., 2010). The two major forms of Arabidopsis glucosinolates are aliphatic glucosinolates (AG) synthesized from methionine, and indolic glucosinolates (IG) synthesized from tryptophan (Sønderby et al., 2010), and their levels can be induced by JA treatment (Kliebenstein et al., 2002). Consistent with the enrichment of the GO term glucosinolate metabolism (GO:0019760) in our omics data ( Figure 1C; Supplemental Data Set 2), the transcript abundance and the levels of the encoding proteins of many IG biosynthetic genes and the genes that are involved in both IG and AG biosynthesis increased in fer relative to WT ( Figure 5, A and B). Consistent with the increased gene expression, we observed a rise in the levels of three different IGs in the fer mutant, indol-3-ylmethyl (I3M), 4-methoxy-indol-3-ylmethyl (4MOI3M), and 1-methoxy-indol-3-ylmethyl (1MOI3M; Figure 5, C-E). Furthermore, while we did not see significant changes in the levels of I3M or 4MOI3M in the nai1 mutant compared to WT, their levels did decrease in the fer nai1 double mutant compared to those measured in fer, suggesting that NAI1 functions downstream of FER to promote the biosynthesis of I3M and 4MOI3M. Disruption of NAI1 did not alter 1MOI3M biosynthesis ( Figure 5E), consistent with previous observations that 1MOI3M is regulated independently from other IGs (Kliebenstein et al., 2002).
To further strengthen the correlation between IG levels and the expression levels of IG genes, we carried out RT-qPCR in 3-week-old plants for six genes involved in IG biosynthesis: CYP79B2, and CYP79B3. Consistent with the increased IG levels in fer, all six genes displayed increased transcript levels in fer ( Figure 6A). While we observed no major changes in the expression of these six genes in the nai1 mutant, the expression of MYC2 decreased in the fer nai1 double mutant, suggesting that NAI1 is required for optimal MYC2 expression and also supports the observation that I3M and 4MOI3M levels are lower in fer nai1 compared to fer. We also performed RT-qPCR using the FERpro:FER-GFP and FERpro:FER K565R -GFP complementation lines. While FERpro:FER-GFP completely restored the expression of all six genes back to WT levels, we witnessed little to partial rescue in FERpro:FER K565R -GFP #5 and partial to complete rescue in FERpro:FER K565R -GFP #8 ( Figure 6B). These results demonstrated that optimal FERK activity is important for the regulation of IG biosynthesis. Taken together, our results suggest that FER negatively regulates indole glucosinolate biosynthesis while inhibiting ER body formation, which strongly corroborates the previous findings of co-regulation between ER body-related genes and glucosinolates biosynthesis and catabolic genes (Nakano et al., 2017; Wang et al., 2017). Further, these results suggest that FER functions as a potential link between environmental signals and ER body formation/IG biosynthesis in response to stress conditions. FER negatively regulates ABA signaling during cotyledon greening through ABI5 FER has been shown to negatively regulate ABA signaling by interacting with ABI2 (Yu et al., 2012;Chen et al., 2016). An interesting finding from our phosphoproteomics analysis was that a group of TFs whose transcription is induced by ABA, including ABF1 (At1g49720), ABF2 (At1g45249), ABF3 (At4g34000), ABF4 (At3g19290), ABA-RESPONSIVE ELEMENT BINDING PROTEIN3 (AREB3; At3g56850), Figure 4 NAI2 functions downstream of FER in ER body formation. A, Immunoblot showing the decreased FER protein levels in the three lines of FER-amiRNA nai2-2 used for the ER body analysis. Ponceau S staining of Rubisco is shown as loading control. B and D, Number of ER bodies in 5day-old (B) and 10-day-old (D) cotyledons from the controls (GFP-HDEL, FER-amiRNA GFP-HDEL #2, and nai2-2 GFP-HDEL) and three individual lines of FER-amiRNA nai2-2 GFP-HDEL (#1, #2, and #3). Data are shown as violin plots with median, first and third quartiles. Different letters indicate significant differences according to one-way ANOVA Tukey's multiple range tests (P 5 0.05). n = 11-35. C and E, Representative confocal images from 5-day-old (C) and 10-day-old (E) cotyledons. Scale bars = 20 mm.
ENHANCED EM LEVEL (EEL; At2g41070), and FBH3 (At1g51140; Song et al., 2016) are hypophosphorylated in fer, suggesting that FER regulates their phosphorylation. Direct target genes were reported for ABF1, ABF2, ABF3, and FBH3 from chromatin immunoprecipitation deep sequencing (ChIP-seq) studies (Song et al., 2016). We observed a significant overlap between fer-regulated genes and the target genes of these ABA-related TFs (Supplemental Figure  S12A; Supplemental Data Set 5), suggesting that FER can function through these TFs to regulate ABA signaling.
Among the peptides that are hypophosphorylated in fer, we mapped the phosphopeptide QG(pS)LTLPR to ABF1, ABF2, ABF3, ABF4, AREB3, and EEL. Further examination revealed that the phosphorylated serine residue and seven out of the eight amino acids in the peptide are also conserved in ABI5, a transcription factor involved in ABAmediated cotyledon greening ( Figure 7A), raising the possibility that FER regulates ABA-mediated cotyledon greening through ABI5. ABI5 can be phosphorylated at many sites by protein kinases such as SNF1-RELATED PROTEIN KINASE2 (SnRK2), BR-INSENSITIVE2 (BIN2), and PROTEIN KINASE SOS2-LIKE5 (PKS5, also named CBL-INTERACTING PROTEIN KINASE11 [CIPK11]), and the kinase(s) responsible for phosphorylating Ser-145 in ABI5 has yet to be identified (Lopez-Molina et al., 2002;Yu et al., 2015). To test the hypothesis that FER regulates ABI5 to regulate cotyledon greening, we generated the fer abi5-7 double mutant. In the cotyledon greening assay, while fer was hypersensitive and abi5-7 was resistant to 1-mM ABA treatment, as quantified by the percentage of seedlings with green cotyledons, the fer abi5-7 double mutant showed a level of tolerance to ABA similar to that of abi5-7, suggesting that FER represses ABI5 function during cotyledon greening (Figure 7, B and C; Supplemental Figure S12B). Further analysis of the cotyledon greening with complementation fer lines showed that FERpro:FER-GFP can largely rescue the ABA hypersensitivity of fer. The complementation by FERpro:FER K565R -GFP was partial, and the degree of the complementation positively correlated with FER K565R -GFP protein abundance, as measured with an anti-GFP antibody (#5 5 #8 5 #6) (Figure 7, D and E; Chakravorty et al., 2018), which demonstrates that FER regulation of ABA-mediated cotyledon greening is dependent on its kinase activity.
To elucidate the potential mechanisms by which FER regulates ABI5, we co-infiltrated the constructs ABI5-FLAG and FER-GFP, mFER-GFP (FER K565R -GFP), or GFP   control into N. benthamiana leaves. Immunoblots established the lower accumulation of ABI5 when ABI5-FLAG was co-infiltrated with FER-GFP, compared to the GFP control ( Figure 7F), suggesting that FER negatively regulates ABI5 protein levels. mFER-GFP did not appreciably change ABI5 protein abundance, suggesting that the negative regulation of ABI5 by FER is kinase dependent. We further generated ABI5-yellow fluorescent protein (YFP) transgenic plants by introducing a transgene consisting of the ABI5 coding sequence cloned in-frame and upstream of YFP and driven by the cauliflower mosaic virus 35S promoter (35Spro:ABI5-YFP OX). We then crossed the resulting transgenic lines to fer and generated fer 35Spro:ABI5-YFP OX lines (Supplemental Figure S12C). To test the effect of FER on ABI5 stability, we treated the rosette leaves of 35Spro:ABI5-YFP OX and fer 35Spro:ABI5-YFP OX with 10 mM ABA. While we detected low levels of ABI5-YFP in ABI5-YFP OX prior to ABA treatment, we observed a modest increase in ABI5-YFP abundance after ABA treatment ( Figure 7G). In contrast, we noticed the greater accumulation of ABI5-YFP in fer ABI5-YFP OX lines even before ABA treatment, which further increased upon ABA treatment. These results suggest that FER negatively regulates ABI5 protein stability ( Figure 7G). An in vitro kinase assay showed that ABI5 can be directly phosphorylated by FER ( Figure 2C). We thus mutated Ser-145 in ABI5 to Ala to produce the nonphosphorylatable variant ABI5 S145A . In vitro kinase assay determined that FER phosphorylation of ABI5 S145A is greatly reduced compared to that of intact ABI5, suggesting that FER phosphorylates ABI5 at S145 ( Figure 7H). Furthermore, we generated transgenic Arabidopsis lines overexpressing ABI5 S145A -YFP FER negatively regulates ABA response during cotyledon greening through ABI5. A, The phosphopeptide QG(pS)LTLPR is shown here, along with a partial alignment between ABI5 and six other ABA-responsive transcription factors in which the phosphopeptide was identified. B, Cotyledon greening assay of WT, fer, abi5-7, and fer abi5-7, on control half-strength MS medium or half-strength MS medium containing 1 mM ABA. Images are of 5-day-old seedlings. C, Quantification of cotyledon greening rates from (B) with 30 seeds of each genotype per treatment. Data represent means ± SEM from three biological replicates. Different letters indicate significant differences according to one-way ANOVA Tukey's multiple range tests (P 5 0.05). D, Cotyledon greening assay of WT, fer, fer FERpro:FER-GFP, or fer FERpro:FER K565R -GFP complementation lines on control half-strength MS medium or half-strength MS medium containing 1 mM ABA. Images are of 5-day-old seedlings. E, Quantification of cotyledon greening rates from (D) with 40 seeds of each genotype per treatment. Data represent means ± SEM from four biological replicates. Different letters indicate significant differences according to one-way ANOVA Tukey's multiple range tests (P 5 0.05). F, Co-infiltration of ABI5-FLAG with FER-GFP or FER K565R (mFER-GFP) in N. benthamiana for 2 days. Proteins were detected by immunoblotting with anti-FLAG (rabbit) and anti-GFP (rabbit) antibodies. Ponceau S staining of Rubisco protein serves as loading control. Data from two independent leaves are shown here. G, Immunoblot showing ABI5 protein abundance upon treatment of ABI5-YFP OX and fer ABI5-YFP OX with 10 mM ABA. Asterisk indicates background bands. Arrow indicates ABI5-YFP. Ponceau S staining of Rubisco protein serves as loading control. H, In vitro kinase assay showing the phosphorylation of ABI5 and decreased phosphorylation of ABI5 S145A and FERK autophosphorylation (top); bottom, SYPRO RUBY staining of the proteins used in the assay. I, Immunoblot showing the stability of ABI5 and ABI5 S145A upon CHX treatment, in transgenic plants overexpressing ABI5-YFP or ABI5 S145A -YFP. Ponceau S staining of Rubisco protein serves as loading control.
(ABI5 S145A -YFP OX). ABI5 S145A -YFP appeared to be more stable than ABI5-YFP when seedlings were treated with the protein synthesis inhibitor cycloheximide (CHX; Figure 7I), suggesting that FER destabilizes ABI5 through phosphorylation at S145 in planta. Similar to ABI5-YFP OX seedlings, ABI5 S145A -YFP OX lines were also hypersensitive to ABA during cotyledon greening (Supplemental Figure S12D). Taken together, our results indicate that FER negatively regulates ABA responses during cotyledon greening through phosphorylation and destabilization of ABI5, in addition to the known regulation of ABI2 by FER (Yu et al., 2012;Chen et al., 2016).

Discussion
The receptor kinase FER plays critical roles in mediating plant growth and development, as well as responses to biotic and abiotic stress. Our integrated omics analysis of the fer-4 mutant not only corroborated previous findings but also revealed new pathways and potential underlying mechanisms of FER function (Figure 8). We showed here that FER negatively regulates ER body formation by regulating the transcription factor NAI1, along with the protein encoded by its target gene NAI2. FER also repressed indole glucosinolate biosynthesis, which was consistent with the previously observed co-occurrence of ER body formation and IG biosynthesis and catabolism (Nakano et al., 2017;Wang et al., 2017). In addition, our phosphoproteomics analysis identified a group of TFs whose encoding genes are induced by ABA and whose phosphorylation is regulated by FER and likely act downstream of FER. We also showed that ABI5 functions downstream of FER and is phosphorylated at the S145 residue by FER, leading to ABI5 destabilization. Thus, our integrated omics study provided new insights into FER functions and underlying molecular mechanisms.
Although ER bodies have been implicated in stress responses, how they are regulated by signaling pathways remains largely unknown. Unlike young seedlings with many ER bodies, the rosette leaves of adult Arabidopsis plants and older seedlings are largely free of ER bodies (Hayashi et al., 2001;Wang et al., 2017). The fact that the fer mutant retains ER bodies in older seedlings and adult leaves indicated that Figure 8 Known and newly identified functions of FER and underlying mechanisms. FER and its co-receptors LLGs/LRE can perceive RALF peptides as ligands (Haruta et al., 2014;Li et al., 2015;Xiao et al., 2019) to regulate diverse biological processes, such as female fertility (Escobar-Restrepo et al., 2007;Hou et al., 2016;Duan et al., 2020), plant growth and development (Guo et al., 2009a;Duan et al., 2010), abiotic and biotic stress responses (Kessler et al., 2010;Chen et al., 2016;Masachis et al., 2016;Stegmann et al., 2017;Feng et al., 2018;Guo et al., 2018;Zhao et al., 2018). Many molecules and FER-interacting proteins are also involved in FER signaling. These include extracellular pectin, reactive oxygen species (ROS), Ca 2 + and nitric oxide (NO), and LRX proteins (Duan et al., 2014;Ngo et al., 2014;Feng et al., 2018;Lin, 2018;Zhao et al., 2018;Dunser et al., 2019;Duan et al., 2020); plasma membrane-associated proteins such as early nodulin-like proteins, Arabidopsis G-protein beta subunit 1, ROPs, FLS2, EFR, and BRI1-associated receptor kinase (Duan et al., 2010;Hou et al., 2016;Stegmann et al., 2017;Yu et al., 2018); intracellular proteins such as ABI2, RIPK, MARIS, MYC2, EBP1, eIF4E1, and intracellular ROS, Ca 2 + , and NO (Yu et al., 2012;Boisson-Dernier et al., 2015;Du et al., 2016;Guo et al., 2018;Li et al., 2018;Zhu et al., 2020). The ligand RALF1 functions through FER to inhibit root elongation (Haruta et al., 2014); FER also interacts with and activates ROP2 to promote root hair development (Duan et al., 2010); FER controls pectin and NO to contribute to female reproduction (Duan et al., 2020); LRXs interact with FER directly or through RALFs to mediate vacuolar expansion or salt tolerance, respectively (Zhao et al., 2018;Dunser et al., 2019); and upon bacterial pathogen attack, FER can serve as a scaffold for FLS2/EFR or phosphorylates and destabilizes MYC2 to promote immune response (Stegmann et al., 2017;Guo et al., 2018). In our current multiomics study, we identified and validated novel functions for FER and underlying molecular mechanisms. FER negatively regulates ER body and indolic glucosinolate biosynthesis through the negative regulation of the transcription factor NAI1, and positively regulates cotyledon greening through the negative regulation of ABI5. FER negatively regulates ER body formation. We further determined that FER inhibits ER body formation through the negative regulation of NAI1, a master transcription factor controlling ER body-related gene expression such as the ER body resident gene NAI2. ER body formation can also be induced by wounding and phytohormones such as JA in Arabidopsis (Hara-Nishimura and Matsushima, 2003;Geem et al., 2019;Stefanik et al., 2020). We showed previously that FER negatively regulates the JA signaling pathway by phosphorylating and destabilizing MYC2, the master transcription factor in the JA pathway (Guo et al., 2018). MYC2 and its homologs MYC3 and MYC4 activate the expression of the ER body gene TRYPTOPHAN SYNTHASE ALPHA CHAIN1, which is involved in JA-induced ER body formation in Arabidopsis adult leaves (Stefanik et al., 2020), which suggests that FER can also regulate ER body formation by regulating the JA signaling pathway. This study thus provides an important link between FER-mediated signaling pathway and ER body formation, and a platform to study the mechanisms by which internal and external stimuli are integrated through the receptor kinase FER to regulate ER body formation and plant stress responses.
Glucosinolates are a group of secondary metabolites that exert important functions in plant responses to biotic stresses (Burrow et al., 2010). ER body-related genes and glucosinolate biosynthetic and catabolic genes are strongly co-expressed, and ER body formation co-occurs with indole glucosinolates (Nakano et al., 2017;Wang et al., 2017). ER bodies are known to accumulate beta-glucosidases such as PYK10 and BGLU21 that are involved in the hydrolysis of glucosinolates and defense against herbivores such as woodlice (Armadillidium vulgare), thereby influencing Arabidopsis/endophyte interactions (Sherameti et al., 2008;Nakazaki et al., 2019;Yamada et al., 2020). In the fer mutant, ER body formation is constitutive with increased levels of beta-glucosidases, as well as elevated indole glucosinolates, which suggests a role for FER in herbivory and other biotic interactions. Future studies should help establish the functions of FER and underlying molecular mechanisms in plant-biotic interactions mediated by ER body and glucosinolates.
Our study further highlights that FER regulates diverse biological processes through the regulation of transcription factors. In addition to our previous finding that FER regulates MYC2 to modulate JA signaling and plant immunity (Guo et al., 2018), here we discovered that (1) FER negatively regulates NAI1 in controlling ER body formation and partially controlling IG biosynthesis, likely through intermediate proteins (Figures 3-6) and (2) that FER phosphorylates and destabilizes ABI5 to control cotyledon greening (Figure 7). In addition, we looked for all 2,492 transcription factors encoded by the Arabidopsis genome (Pruneda-Paz et al., 2014) among differentially expressed transcripts, DAPs, and phosphoproteins in fer. This analysis showed that more than 20% of the transcription factors (524/2,492) exhibit altered levels for their transcripts, encoded proteins, or phosphorylation state (Supplemental Figure S13; Supplemental Data Set 6), which strongly suggests that FER regulates diverse biological processes through the extensive involvement of multiple transcription factors, directly or indirectly. Future analysis of these transcription factors should reveal a more comprehensive FER signaling network and provide a better understanding of the molecular interplay in FER-mediated signaling pathways.

Protein extraction and digestion
The proteomics experiments were carried out based on established methods as follows Song et al., 2018aSong et al., , 2018bSong et al., , 2020. Three biological replicate samples were collected from 4.5-week-old entire rosettes from both Col-0 and the fer mutant, with 4-5 Col-0 rosettes and 8-10 fer rosettes per replicate. Lysis buffer consisting of 8 M urea, 100-mM Tris-HCl pH 7, 5-mM tris(2-carboxyethyl)phosphine (TCEP) and 1 Â phosphatase inhibitor (2.5-mM NaF, 0.25-mM NaVO 4 , 0.25-mM sodium pyrophosphate decahydrate, and 0.25-mM glycerophosphate in H 2 O) was added to 250mg tissue at a ratio of 1:2 sample:buffer (w:v). Zirconium oxide beads (1-mm diameter, Next Advance) were added to the samples at a 1:1 ratio (v:v) and then the samples were shaken using a GenoGrinder (SPEX) at 1,500 rpm for 3 min. The samples were centrifuged at 4,000 g for 3 min. The shaking and centrifugation steps were repeated once. Samples were transferred to new tubes, to which four volumes of prechilled 100% acetone were added. Samples were precipitated at -20 C for 430 min followed by centrifugation at 4,500 g for 10 min at 4 C. Acetone (80%, v/v) was added to the pellets and the samples were probe-sonicated to resuspend the pellet and shear DNA. Samples were incubated at -20 C for 45 min and then centrifuged at 4,500 g for 10 min at 4 C. Precipitation and sonication in 80% (v/v) acetone was performed 3 times in total. Prechilled 100% methanol was then added to each pellet, the samples were probe-sonicated, and kept at -20 C for 30 min prior to centrifugation at 4,500 g for 10 min at 4 C. Methanol precipitation was repeated once. The supernatant was discarded and the pellet was placed in a vacuum concentrator until nearly dry.
Proteins were solubilized in 0.5-mL protein resuspension buffer (8 M urea, 0.1-M Tris-HCl pH 7.0, 5-mM TCEP, 1 Â phosphatase inhibitor) and probe-sonicated. The protein amount was evaluated by Bradford assay and $1 mg was mixed with 3.5-mL urea solution (8 M urea, 0.1-M Tris-HCl pH 8.0, 1 Â phosphatase inhibitor). This solution was added to an Amicon Ultracel-30K centrifugal filter (Cat # UFC803008) and centrifuged at 4,000 g for 20-40 min at room temperature. This step was repeated once. Then 4 mL of urea solution with 2-mM TCEP was added to the filter unit and centrifuged at 4,000 g for 20-40 min. Next, 2-mL iodoacetamide (IAM) solution (50-mM IAM in 8-M urea) was added and incubated without mixing at room temperature for 30 min in the dark prior to centrifugation at 4,000 g for 20-40 min. Two milliliters of 8-M urea were added to the filter unit, which was then centrifuged at 4,000 g for 20-40 min. This step was repeated once. Then, 2 mL of 0.05 M NH 4 HCO 3 with 1 Â phosphatase inhibitor was added to the filter unit and centrifuged at 4,000 g for 20-40 min. This step was repeated once. Then the filter unit was added to a new collection tube and 2 mL of 0.05 M NH 4 HCO 3 with 1 Â phosphatase inhibitor with trypsin (enzyme to protein ratio 1:100, w/w) was added. Samples were incubated at 37 C overnight. Undigested proteins were estimated using Bradford assays before adding trypsin (1 lg/lL) at a ratio of 1:100 (w/w) and an equal volume of Lys-C (0.1 lg/lL) was added and incubated for an additional 4 h at 37 C. The filter unit was centrifuged at 4,000 g for 20-40 min. One milliliter of 0.05-M NH 4 HCO 3 with 1 Â phosphatase inhibitor was added and centrifuged at 4,000 g for 20-40 min. The samples were acidified to pH 2-3 with 100% formic acid and centrifuged at 21,000 g for 20 min. Finally, samples were desalted using 50 mg Sep-Pak C18 cartridges (Waters, Milford, MA, USA). Eluted peptides were dried using a vacuum centrifuge (Thermo Fisher Scientific, Waltham, MA, USA) and resuspended in 0.1% (v/v) formic acid. Peptide amount was quantified using the Pierce BCA Protein assay kit.

TMT labeling
TMTsixplex label reagents (Thermo Fisher Scientific, Lot #SH254566) were used to label the samples according to the manufacturer's recommended peptide-to-tandem mass tag (TMT) reagent ratio. Four hundred micrograms of vacuumdried peptides from each sample were resuspended with 400-mL 50 mM TEAB buffer and vortexed for 10 min at room temperature. Then, 41-mL acetonitrile was added to each tube of TMT label (0.8 mg), vortexed, and incubated at room temperature for 5 min to resuspend the labels. Four tubes (4 Â 41 mL) of each type of TMT label (3.2-mg TMT label in total) were added to each tube of peptides (400 mg), pipetted up and down several times, and vortexed to mix them well. After a 2-h incubation at room temperature, 32 mL of 5% (w/v) hydroxylamine was added to each tube, vortexed, and incubated at room temperature for 15 min to quench the labeling reaction. Next, the six samples were mixed together, a 35-mg aliquot of peptides was reserved for protein abundance profiling, and the remaining peptides were used for phosphopeptide enrichment, and stored at -80 C.

Phosphopeptide enrichment
The TMT-labeled phosphopeptides were enriched using Titansphere Phos-TiO 2 beads (GL Sciences 5010-21315) based on previously published methods (Kettenbach and Gerber, 2011;Song et al., 2018b). The beads were prepared by resuspending in 1.5-mL wash and binding buffer (2 M lactic acid in 50% [v/v] acetonitrile), vortexing, and then centrifugation at 3,000 g for 1 min at room temperature; this was repeated a total of 3 times. At the last washing step, 5 and 11 mg of TiO 2 beads were aliquoted into new tubes before centrifugation. After centrifugation, the wash and binding buffer were removed and the TiO 2 beads were saved for phosphopeptide enrichment. About 2.4 mg of TMT6-labeled and vacuum-dried peptides were resuspended with 2.4 mL wash and binding buffer and then added to the tube containing 11 mg TiO 2 beads, rotated at room temperature for 1 h, and then centrifuged at 3,000 g for 1 min at room temperature. The supernatant was processed with a second round of enrichment using 5 mg of TiO 2 beads. Then, 1.8-mL wash and binding buffer was added to each tube from the two enrichment steps, vortexed, and centrifuged at 3,000 g for 1 min at room temperature. This wash was repeated once. Next, the TiO 2 beads were washed twice with 1.8 mL of 50% (v/v) acetonitrile in 0.1% (w/v) trifluoroacetic acid. After the wash steps, 500 mL of 3% (w/v) ammonium hydroxide was added to each tube of the two enrichment steps, vortexed and centrifuged at 3,000 g for 1 min at room temperature. The eluted supernatants were combined. One more elution step was performed with 5% (w/v) ammonium hydroxide. All supernatants from the two elution steps were combined and dried in a speed vac; the phosphopeptides were then resuspended in 0.1% (w/v) FA, and stored at -80 C until liquid chromatography-tandem mass spectrometry (LC/MS-MS) run.
Eluted peptides were analyzed using a Thermo Scientific Q-Exactive Plus high-resolution quadrupole Orbitrap mass spectrometer, which was directly coupled to the HPLC. Data-dependent acquisition was obtained using Xcalibur version 4.0 software in positive ion mode with a spray voltage of 2.00 kV and a capillary temperature of 275 C and a retention factor of 60. MS1 spectra were measured at a resolution of 70,000, an automatic gain control (AGC) of 3 Â 10 6 with a maximum ion time of 100 ms and a mass range of 400-2,000 m/z. Up to 15 MS2 were triggered at a resolution of 17,500 with a fixed first mass of 120 m/z for the phosphoproteome and 115 m/z for the proteome. An AGC of 1 Â 105 with a maximum ion time of 50 ms, an isolation window of 1.3 m/z for phosphoproteome and 1.2 m/z for proteome, and a normalized collision energy of 31 and 32 were used for nonmodified and phosphoproteomes, respectively. Charge exclusion was set to unassigned, 1, 5-8, and 48. MS1 that triggered MS2 scans were dynamically excluded for 25 or 30 s for nonmodified and phosphoproteomes, respectively.

Data analysis
The raw data were analyzed using MaxQuant version 1.6.1.0 (Tyanova et al., 2016). Spectra were searched, using the Andromeda search engine (Cox et al., 2011) against the Arabidopsis TAIR10 proteome file entitled "TAIR10_ pep_20101214" that was downloaded from the TAIR website (https://www.arabidopsis.org/download/index-auto.jsp? dir=%2Fdownload_files%2FProteins%2FTAIR10_protein_lists) and was complemented with reverse decoy sequences and common contaminants by MaxQuant. Carbamidomethyl cysteine was set as a fixed modification, while methionine oxidation and protein N-terminal acetylation were set as variable modifications. For the phosphoproteome, "Phosho STY" was also set as a variable modification. The sample type was set to "Reporter Ion MS2" with "6plex" "MT" selected for both lysine and N termini. TMT batch-specific correction factors were configured in the MaxQuant modifications tab (TMT Lot SH254566). Digestion parameters were set to "specific" and "Trypsin/P;LysC." Up to two missed cleavages were allowed. A false-discovery rate, calculated in MaxQuant using a target-decoy strategy (Elias and Gygi, 2007) of less than 0.01 at both the peptide spectral match and protein identification level was required. The "second peptide" option to identify co-fragmented peptides was not used. The match between runs feature of MaxQuant was not utilized.
Prior to statistical analysis, protein/phosphosite intensity data were normalized such that the total intensity for each TMT lane was equal across the run (referred to as sample loading normalization). No imputation for missing values was performed. Then, statistical analysis was performed using the PoissonSeq package in R . Proteins were categorized as differentially abundant if they had a q 4 0.05. For phosphosites, we categorized differential abundance as those sites with q 4 0.2 based on the q-value histogram (Supplemental Figure S1).

RNA extraction and 3 0 -RNA sequencing
Total RNA was extracted using RNeasy Plant Mini Kit (Qiagen, Hilden, Germany; catalog # 74904), and genomic DNA contamination was removed using RNase-free DNase I Set (Qiagen; catalog # 79254) in-column during RNA extraction according to the manufacturer's protocols. Libraries were constructed using QuantSeq 3 0 -mRNA-Seq Library Prep Kit from Illumina. Libraries were run on a HiSeq4000 with 50-bp reads, as QuantSeq is optimized for 50-to 100-bp reads (Moll et al., 2014). Each library was run twice to increase the average read count to $6 million reads per sample. Read alignment to the TAIR10 genome was performed using the STAR aligner (Dobin et al., 2013) using default parameters exce--outFilterMismatchNoverLmax was set to 0.6 to account for the shorter QuantSeq reads, resulting in $70% of reads uniquely mapped per sample, which is expected for QuantSeq (Moll et al., 2014). Counts per transcript were obtained using htseq-count using the intersection-nonempty parameter for reads spanning more than one feature (Anders et al., 2015).
GO enrichment analysis and network reconstruction were performed using the ClueGO application in Cytoscape (Bindea et al., 2009). Terms were considered enriched if they had a corrected P 50.05.
The enrichment of consensus sequences for FER phosphorylation was performed using the motifeR R package (Wang et al., 2019a(Wang et al., , 2019b with default settings. Sequence logos were made using the ggseqlogo R package (Wagih, 2017). Both analyses were performed in R version 3.6.2 (R Core Team, 2019).
The sequences encoding the substrate proteins: ABF3, OXS2, ABI5, ABI5 S145A , SHOU4, and NAIP1 were PCR amplified and cloned into pMBP-H to generate MBP fusion proteins (Guo et al., 2018). The cloning primers are listed in Supplemental Table S1. The coding sequences of FBH3, GT2, NSI, and At35g1950 cloned into pDEST22 were obtained from TAIR (https://www.arabidopsis.org/). The coding sequences were then cloned into pDONR-221 and then pDEST-HIS MBP to generate MBP fusion proteins. All the constructs used in this study were listed in Supplemental  Table S2.

In vitro kinase assay
The in vitro kinase assay was carried out as described (Guo et al., 2018). Briefly, recombinant purified GST-FERK was mixed with MBP-substrate proteins in 20-mL reactions, with 10 mCi 32 P-cATP and incubated at room temperature for 1 h. The reactions were stopped using 20 mL 2 Â SDS sample buffer and resolved on an 8% (w/v) SDS-PAGE.

BiFC assays
The mFER coding sequence was cloned in-frame and upstream of the sequence encoding the N terminus of YFP, while the SHOU4, FBH3, ABI5, and DDL coding sequences were cloned in-frame and upstream of the sequence encoding the C terminus of YFP. All resulting plasmids were individually introduced into Agrobacterium (Agrobacterium tumefaciens) (strain GV3101). Positive Agrobacterium colonies were cultured in liquid LB medium containing 0.2-mM acetosyringone for 1-2 days. Cells were collected by centrifugation, washed, and resuspended in infiltration buffer (10-mM MgCl 2 , 10-mM MES, pH 5.7, 0.2-mM acetosyringone) to a final optical density (OD) measured at a wavelength of 600 nm (OD 600 ) of 0.9. Agrobacteria carrying appropriate pairs of nYFP and cYFP constructs were mixed in a 1:1 ratio and infiltrated into the lower surface of N. benthamiana leaves from 2-month-old plants grown on soil. After 36 h, YFP signals were detected using a Zeiss Laser Scanning Microscope 700 (LSM700).

Transient infiltration assays
Agrobacteria carrying the constructs of interest were grown in liquid LB medium with antibiotics in a 30 C shaker for 2 days. After collecting the agrobacteria by centrifugation, the cells were resuspended in infiltration buffer (10-mM MgCl 2 , 10-mM MES pH 5.7, 200-lM acetosyringone) to a final OD 600 of 0.3. Leaf infiltration was conducted with a 1-mL syringe on the abaxial side of the leaves, as above. At least two biological replicates were examined for each target construct.
Two days after infiltration, five leaf discs (7 mm in diameter) were collected for each sample and flash-frozen in liquid nitrogen and ground directly in 200 lL of 2 Â SDS sample buffer. The samples were resolved on 8% (w/v) SDS-PAGE, followed by immunoblotting using commercial rabbit anti-GFP (A11122; Invitrogen, Waltham, MA, USA) or anti-FLAG antibodies (F7425; Sigma-Aldrich, St Louis, MO, USA) at a 1:1,000 dilution.

Co-IP
Agrobacteria carrying constructs encoding FLAG-and GFPtagged proteins of interest were co-infiltrated into N. benthamiana leaves. Leaf samples were collected 2 days after the infiltration. One gram of each sample was ground in liquid nitrogen and homogenized in 2.5 mL IP buffer (10-mM HEPES Ph 7.5, 100-mM NaCl, 1-mM EDTA, 10% [v/v] glycerol, and 0.5% [v/v] Nonidet P-40) with 1-mM PMSF, 20-mM MG132, and proteinase inhibitor cocktail for 20 min at 4 C with rotation. Five micrograms of FLAG M2 antibody (F1804; Sigma-Aldrich) was prebound to 60-mL protein G Dynabeads (10003D; Thermo Fisher Scientific) for 30 min in 1Â phosphate buffer saline (PBS) buffer with 0.02% (v/v) Tween-20 at room temperature. The beads were washed once with the same PBS buffer and resuspended in Co-IP buffer. After protein extraction, 20 mL of anti-FLAG prebound Dynabeads was added to each sample for another 1.5 h incubation at 4 C with rotation. Dynabeads were precipitated using DynaMagnetic rack (12321D; Thermo Fisher Scientific) and washed twice with Co-IP buffer with 0.5% (v/v) Nonidet P-40 and twice with Co-IP buffer without Nonidet P-40. The IP product was resuspended in 2Â SDS sample buffer and used for immunoblotting with rabbit anti-GFP (A11122; Invitrogen) and rabbit anti-FLAG antibody (F7425; Sigma-Aldrich).

ER body observation and quantification
For ER body observation, GFP-HDEL-labeled structures were observed and images were taken by confocal microscopy using a Zeiss Laser Scanning Microscope 700 (LSM700) with a 63Â oil immersion objective. GFP was excited with a 488 nm laser line and detected at 555 nm. The ER body quantification was carried out as described (Yamada et al., 2008).

Measurement of glucosinolate contents
Glucosinolates were measured as previously described (Brown et al., 2003;Kliebenstein et al., 2001a-c). Briefly, 3week-old whole rosettes were pooled, weighed, frozen, and harvested in 90% (v/v) methanol. Tissues were homogenized for 3 min in a paint shaker, centrifuged, and the supernatants were transferred to a 96-well filter plate with DEAE Sephadex. The filter plate with DEAE Sephadex was washed once with water, once with 90% (v/v) methanol, and once with water again. The Sephadex-bound glucosinolates were eluted after an overnight incubation with 110 lL of sulfatase. Individual desulfo-glucosinolates within each sample were separated and detected by high-performance liquid chromatography (HPLC)-diode array detection, identified using retention time and absorbance spectra developed from purified compounds and re-validated using Arabidopsis genotypes with known chemotypes, quantified using relative response factors developed by comparison to standard curves from purified compounds as previously reported (Brown et al., 2003), and further normalized to fresh weight.

Cotyledon greening assay
Surface-sterilized seeds were germinated on control halfstrength MS medium or half-strength MS medium containing 1-mM ABA under constant light. Cotyledon greening was observed 4-5 days later. The cotyledon greening rate was calculated using the percentage of seeds with green cotyledons out of all seeds sown, for each genotype.

CHX and ABA treatments for protein stability assays
For ABI5 protein stability assays, 3-week-old rosette leaves were collected and cut into smaller pieces and incubated in liquid half-strength MS medium for 2 h before adding 10-mM ABA. Samples were then collected at the indicated times. Time 0 was collected right before the addition of ABA. Total proteins were extracted using 2Â SDS buffer and samples were immunoblotted with rabbit anti-GFP antibody (A11122; Invitrogen).
For ABI S145A -YFP protein stability assays, 7-day-old seedlings of the ABI5-YFP and ABI S145A -YFP transgenic lines were treated with 0.5-mM CHX in half-strength liquid MS medium for the indicated times. Total proteins were extracted in 2Â SDS loading buffer and samples were immunoblotted with rabbit anti-GFP antibody (A11122; Invitrogen).
For endogenous FER protein detection, samples were immunoblotted with lab-made rabbit anti-FER antibody.

RT-qPCR
Total RNA was isolated from 3-week-old rosette leaves using an RNeasy Plant Mini Kit (Qiagen; catalog # 74904) and genomic DNA was removed using RNase-free DNase Set (Qiagen; catalog # 79254) in-column during RNA extraction. The first-strand cDNA was synthesized with an iScript cDNA Synthesis Kit (BioRad, Hercules, CA, USA; catalog # 1708891). Real-time PCR was performed using SYBR Green PCR Master Mix (Applied Biosystems, Waltham, MA, USA; catalog # 4309155) on the StepOnePlus Real-Time PCR system (Applied Biosystems). Relative gene expression was determined by applying the 2 -DDCT (CT, cycle threshold) method and normalized to the expression of the reference gene ACTIN2 (At3g18780). RT-qPCR was performed with two to three technical replicates from three to four independent biological replicates. Primers used for qPCR are provided in Supplemental Table S1.

Statistical analysis
Graphs were created in GraphPad Prism software (version 9.3.0). SPSS version 27.0 software (IBM, Armonk, NY, USA) was used for statistical data analysis. The data are shown as means ± standard error of the mean (SEM) and were subjected to one-way analysis of variance (ANOVA) Tukey's multiple range tests (P 5 0.05). ANOVA data are provided in Supplemental Data Set 7.

Supplemental data
The following materials are available in the online version of this article.
Supplemental Figure S2. The comparisons of differentially expressed transcripts and proteins in fer, and the differentially expressed transcripts of this study to that of the previous publication.
Supplemental Figure S3. Enriched GO terms in transcripts with increased levels in fer.
Supplemental Figure S4. Enriched GO terms in transcripts with decreased levels in fer. Figure S5. Enriched GO terms in proteins with increased levels in fer.

Supplemental
Supplemental Figure S6. Enriched GO terms in proteins with decreased levels in fer.
Supplemental Figure S7. Enriched GO terms in phosphoproteins with increased levels in fer.
Supplemental Figure S8. Enriched GO terms in phosphoproteins with decreased levels in fer.
Supplemental Figure S9. FER directly phosphorylates many proteins with diverse functions.
Supplemental Figure S11. FER negative regulation of ERbody formation is validated using FER-miRNA knockdown.
Supplemental Figure S13. FER regulates 420% of TFs in the Arabidopsis genome, directly or indirectly.
Supplemental Table S1. Primers used in this study. Supplemental Table S2. Constructs used in this study. Supplemental Data Set 1. Misregulated gene list in FER omics data.
Supplemental Data Set 2. Enriched GO terms in FER omics data.
Supplemental Data Set 3. Enriched consensus sequences of FER phosphorylation sites.
Supplemental Data Set 6. Transcription factors that are differentially expressed in FER omics data.
Supplemental Data Set 7. ANOVA results in this study.